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Abstract 

Lean, prevaporized, premixed combustors 
are susceptible to combustion-acoustic 
instabilities. A model was developed to predict 
eigenvalues of axial modes for combustion- 
acoustic interactions in a premixed combustor. 
This work extends previous work by including 
variable area and detailed chemical kinetics 
mechanisms, using the code LSENS. Thus the 
acoustic equations could be integrated through 
the flame zone. Linear perturbations were made 
of the continuity, momentum, energy, chemical 
species, and state equations. The qualitative 
accuracy of our approach was checked by 
examining its predictions for various unsteady 
heat release rate models. Perturbations in fuel 
flow rate are currently being added to the model. 


Introduction 

Lean, prevaporized, premixed combustors 
are being considered for future low NO x engines. 
These combustors are susceptible to instabilities 
due to feedback between pressure fluctuations 
and combustion, causing damaging mechanical 
vibrations of the system, as well as degrading 
emissions characteristics and combustion 
efficiency. In a lean flame, blowout can also 
occur. 

Many studies have been made on the 
interactions between acoustics and combustion 
for various applications, including liquid rockets, 
solid rockets, industrial gas turbines, aircraft gas 
turbines, pulse combustors, afterburners, and 
ramjet combustors; see Oyediran et al\ The 
current analysis adds to previous work by 
including detailed reaction mechanisms, multiple 


reflections and transmissions at the flame, and 
variable area. 

The acoustics of the combustion chamber 
become a problem when the coupling of 
fluctuating heat release rate and pressure 
perturbations is concomitant with insufficient 
damping in the chamber and transmission at the 
boundaries to eliminate the amplification of 
acoustic energy due to the flame. Many factors 
affect the fluctuating heat release rate: (1) 
turbulence in the incoming flow, (2) fluctuations in 
the fuel delivery system (either caused by its own 
dynamics, or coupled with pressure fluctuations 
in the chamber), (3) variations in atomization and 
vaporization of the fuel and in fuel-air mixing, (4) 
hydrodynamic instability and vortex shedding at 
the flame holder, (5) variations in chemical kinetic 
reaction rates, (6) mean heat release rate, and 
(7) combustor type and geometry (e.g., premixed 
dump combustors and airblast spray 
combustors). All these factors are potentially 
important in a dynamic heat release rate model. 

The turbulence in the incoming flow is 
invariably amplified by the flame (see Giammar 
and Putnam 2,3 ). The noise tends to be broad 
band and does not necessarily cause acoustic 
resonances, although it could be the source of 
the oscillations amplified when an acoustic 
resonance does develop. 

Fluctuations in the fuel flow rate into the 
flame have a direct effect on the fluctuating heat 
release rate. The fuel injection system may have 
its own resonant frequency, acting as a driving 
function to the heat release rate, or it may couple 
with the dynamics of the combustion chamber, in 
which case the fuel injection system is often 
represented by an admittance. The dynamics of 
the fuel delivery system is a concern for stability 
analyses of combustors 4 . 
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In a combustor with a liquid reactant, the 
instantaneous rate at which reactants are 
delivered to the flame also depends upon the 
atomization, vaporization, and mixing rates. For 
liquid rockets, especially, these rates are 
important, because they often limit the burning 
rate 56 . 

The dynamics of vortices and recirculation 
zones can affect the heat release rate, and are 
particularly important for premixed dump 
combustors. Large variations in the recirculation 
zones can cause separation of the boundary 
layer in the premixing region, leading to 
flashback 7 . 

The heat release rate can also change due 
to variations in the chemical kinetics of the 
combustion reaction. This factor becomes 
important when the time scale of the chemistry is 
of the same order as that of the chamber 
acoustics, or if the heat release is distributed over 
a length scale similar to that of the acoustics. 

Which of the above effects become 
important is determined by the distribution of 
mean heat release rate, combustor type, and fuel 
type (i.e., gaseous or liquid). 

The boundary conditions of the system are 
also important, because they determine how 
much of the acoustic energy leaves the 
combustion chamber. For rocket systems the 
admittance of a choked nozzle can be 
calculated 8 . However, for other systems, the 
acoustic boundary conditions are often not so 
well defined. This is an area requiring additional 
investigation for gas turbine combustors. The 
present work, however, concentrates on the 
dynamics inside the combustion chamber and 
does not address boundary conditions in detail. 

Model 

A linear model was developed to predict the 
dynamics of a one-dimensional chemically 
reacting flow. The model can predict eigenvalues 
of acoustic resonances and standing wave 
patterns for resonance conditions. The 
advantages of the model include the use of 
detailed chemical reaction mechanisms and the 
ability to integrate the acoustic equations through 
the flame zone. Thus the needs to assume a thin 
reaction zone and specify phase and amplitude 
shifts at the flame front are eliminated. Instead, 


these quantities are simply calculated while 
integrating through the flame. 

Solution variables are decomposed into 
mean (steady) and fluctuating components, 


p = p + p', T = T + T', Q = Q + Q', etc. 

(1) 

where the overbar and prime denote mean and 
fluctuating quantities, respectively, and p is the 
pressure, T the temperature, and Q the heat 
release rate per unit length. 

Following the approach of Bloxsidge et al. 9 , 
the one-dimensional mean flow equations are as 
follows: 

Continuity Equation 

-iL(p uA) = 0. (2) 

dx 

Momentum equation 

-~6v + dp = 0 ( 3 ) 

dx dx 

Energy Equation 



In these equations x is the axial location, p the 
density, u the velocity, A the cross-sectional area, 
and c p the constant-pressure, mass-specific heat 
of the mixture. 

The fluctuating components were further 
separated into (1) a function of time t (e**, where 
i = PT and co is the complex oscillation rate), 
to represent the frequency and growth/decay rate 
of oscillations and (2) a function of position, p(x). 
The latter function was itself decomposed into 
real and imaginary parts, to represent the phase 
and magnitude of the oscillations at each spatial 
location, for example, 

p' = pMe 10 ”, 

where 

p(x) = p,(x) + ip 2 (x). (6) 
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The functions p,(x) and p 2 (x) were determined by 
integrating the perturbed equations of state and 
mass, momentum, and energy conservation, 
equations 8-10. 

The temperature, pressure, density, 
velocity, and heat release rate were all perturbed, 
as described above. The perturbation in heat 
release rate, Q', is the energy source for 
amplification of acoustic waves as they propagate 
through the combustion zone: 

Q = Q + Q' = Q + Qe w . ( 7 ) 

The unsteady equations of mass, momentum, 
and energy conservation for linear perturbations 
of frequency co are as follows: 

j!_(puA + puA) = -icopA, (8) 


(pu 


du 

+ pu)_ 
dx 


— du 

+ P u -^“ 

dx 


dp . 

"dx = Hcopu ’(9) 


and 


(pu + pu)[c p T + lu 2 ^ 


+ pu(c p T + uu) 
Q - icoA^pfc v T + 2.U 2 + p(c v T + uu) . 


(10) 


The heat release rate was modeled as 
follows. At a given axial location the heat release 
rate per unit length is a function of the 
temperature, density, and composition at that 
point: 


0=0+0' * f(T, p, a,), i = 1.....NS 


where o, is the mass-specific mole number of 
species i, that is, number of moles of species i 
per unit mass of mixture, and NS is the number 
of chemical species (reacting and inert). The 
perturbation in heat release rate was given by 


Q 


^3Q^ 


/ _ 


3T 


T' + 


/p 


3Q 


K dp jt 


(12) 


3Q/3p are functions of position only, and not time. 
Therefore, they need to be calculated only once 
at each axial location for a given mean flow. 
These partial derivatives were determined using 
a modified version of LSENS, the Lewis Kinetics 
and Sensitivity Analysis Code 10 ". 

In addition to solving the ordinary 
differential equations (ODE's) for the mean 
reacting flow, LSENS generates the linear 
sensitivity coefficients {3y/3T|j} and {3y/3tij}, where 
y, is a dependent variable (species mole number, 
temperature, etc.), y, = dy/dt, and rij is an 
independent parameter of interest (e.g., a rate 
coefficient parameter or an initial condition value). 
The sensitivity coefficients are obtained by 
solving their governing ODE's, by using the 
backward differentiation formula method, as 
implemented in the packaged code LSODE 12 . Of 
particular relevance to the present work are the 
sensitivity coefficients {dq/3T|j}, where q is the 
heat release rate per unit mass: 

NRS 

4- -EM- (13) 

i-1 


Here NRS is the number of reacting species and 
h-, the molar-specific enthalpy of species i. In the 
present model t|, = T and ti 2 = p. If the effects of 
perturbations in fuel injection rate are also 
considered, equation (12) must include sensitivity 
coefficients with respect to the upstream fuel 
concentration, since it will change the local 
composition. The partial derivative 3q/3r|j, 
obtained simply by differentiating equation (13) 
with respect to ify is given by 


-. . NRS 

is. - -t £% 

drii w drii 


NRS 


v-' • 3T 
X>i c P.i=j-> 

i-1 oT|j 


(14) 


where c pi is the constant-pressure, molar-specific 
heat of species i. Equation (14) is used to 
calculate the partial derivatives in equation (12), 
which, in turn, gives the heat release rate 
perturbation in equation (10). 

Results 


The linear stability model described in the 
previous section was used to determine complex 
eigenvalues for a simple combustor. These 
calculations were made in order to test the 


A benefit of this approach is that for linear 
perturbations, the partial derivatives 3Q/3T and 
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qualitative agreement of the model with 
combustion-acoustic theory. In our formulation, 
the real part of the eigenvalue represents the 
oscillation rate and the imaginary part the 
growth/decay rate. A positive imaginary 
component indicates disturbances will decay (i.e., 
a stable system), while a negative imaginary 
component indicates disturbances will grow (i.e., 
an unstable system). 

We first present results for a constant-area 
combustor, and then examine the effects of 
variable area. All results presented herein were 
obtained with a premixed propane-air system and 
a fuel-air equivalence ratio of 0.9. The propane- 
oxygen combustion chemistry was modeled with 
the mechanism of Jachimowski 13 . The reaction 
mechanism consisted of 85 reactions among 28 
species, including the inert species N 2 (i.e., 
reactions involving N 2 were not considered). 

The inlet state and flow geometry were as 
follows. The temperature and pressure were 
1000 K and 5 atm., respectively. The combustor 
length was 2 m and the inlet velocity 100 m/s. 
The flame was stabilized by assuming that a 
flameholder would increase the residence time, 
thereby enabling the combustion reaction to 
occur. However, the full area (i.e., ignoring 
flameholder blockage) and resultant velocity were 
used for the acoustic calculations. The flow 
conditions used in this work are summarized in 
Table 1. 

The present study used the boundary 
conditions u' = 0 and dp'/dx = 0 at the walls. If 
the mean flow Mach number were small, these 
boundary conditions would represent acoustically 
hard walls at the boundaries. For a nonnegligible 
Mach number, however, dp'/dx = 0 is not the 
condition for complete reflection, for which the 


TABLE 1.— CONDITIONS FOR ONE- 
DIMENSIONAL FLOW COMBUSTION- 
ACOUSTIC CALCULATIONS 


Fuel 

Propane 

Flame Type 

Premixed 

Equivalence ratio 

0.9 

Inlet temperature 

1000 K 

Inlet pressure 

5 atm. 

Inlet velocity 

100 m/s 

Combustor length 

2 m 


magnitude of dp'/dx of the wave traveling with the 
flow will be smaller than that of the wave 
traveling against it (for the same amplitude and 
phase angle). 

The model can accept any user-supplied 
admittance as a boundary condition. However, 
the goal of our work was to show that the model 
is in qualitative agreement with theory. Therefore, 
the boundary conditions are not particularly 
important for the present study. They will, 
however, become important when the model 
predictions are compared with experimental data. 

The grid spacing (i.e., step size) used to 
integrate the ODE's is automatically selected by 
LSODE. The variation of the spacing with axial 
position is shown in Figure 1 for the conditions 
given in table 1 and constant area. Note that the 
grid spacing was small in regions where the 
mean flow quantities varied rapidly. The step size 
pattern produces a chemical kinetics solution that 
satisfies prescribed accuracy requirements 10 - 12 , 
thereby enabling the linear perturbation model to 
integrate the governing ODE's through the flame 
region. Figures 2 and 3 show the mean flow 
temperature and Mach number distributions, 
respectively, for the conditions given in table 1 . 

Various models were used for the 
fluctuating heat release rate, Q', in order to 
validate the results qualitatively. It was first set 
equal to zero, to determine the eigenvalues for 
the passive system (i.e., the flame does not add 
acoustic energy). Then, the fluctuating heat 
release rate was modeled (1) as a real constant 



Axial Position [m] 

Figure 1 . — Grid spacing selected by chemical kinetics code 
LSENS for conditions given in table 1 and constant 
area. 
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Axial Position [m] 


Figure 2— Axial distribution of mean flow temperature for 
conditions given in table 1 and constant area. 



Axial Position [m] 


Figure 3.— Axial distribution of mean flow Mach number for 
conditions given in table 1 and constant area. 


times the pressure perturbation (i.e., in phase 
with pressure), (2) as an imaginary constant 
times the pressure perturbation (i.e., n/2 radians 
out of phase with pressure), (3) as a real 
constant times the velocity perturbation (i.e., in 
phase with velocity), and (4) using the 
perturbations in the chemistry from equation 1 2, 
as described in the previous section. 

Eigenvalues were calculated for these heat 
release rate models; their real and imaginary 
components are shown in Figure 4. For each 
point in this figure, table 2 gives the 
corresponding fluctuating heat release rate 
model; here c, and c 2 are real positive constants. 

Points A and B represent the two lowest 
frequency axial modes for the passive system 
using the mean flow conditions shown in Figures 
2 and 3. The standing wave pattern for pressure 
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Figure 4. — Eigenvalues for flow conditions given in table 1 , 
constant area, and dynamic heat release rate models 
described in table 2. 


TABLE 2.— FLUCTUATING HEAT RELEASE 
RATE MODELS USED IN FIGURES 4 AND 8 


Point 

Fluctuating Heat 
Release Rate Model 

A 

Q7Q = 0 

B 

Q7Q = 0 

C 

Q7Q = Ct p' 

D 

Q'/Q = -j c, p' 

E 

Q7Q = i c, p' 

F 

Q7Q = c 2 u' 

G 

Q7G — > Eq. 12 

H 

Q7Q -4 Eq. 12 

1 

Q7Q = 0 

J 

Q7Q = 0 

K 

Q7Q -4 Eq. 12 

L 

Q7Q -> Eq. 12 


and velocity for the first mode (point A) are 
shown in Figure 5. This mode is nominally a half 
wavelength mode. However, since the mean flow 
has a significant velocity relative to the speed of 
sound and since there is a phase shift in the 
standing velocity wave near the flame, the 
standing wave is not exactly a half wavelength of 
the traveling waves. 

The imaginary part of to represents the 
growth/decay rate of the resonance. For the 
passive case (point A) the imaginary part of w is 
210 s' 1 , indicating that the oscillation is decaying. 
This is an important test of the model, since a 
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Axial Position [m] 

Solid Line: normalized pressure 
Dashed Line: normalized velocity 

Figure 5. — Standing wave pattern for lowest frequency mode 
for constant area and passive case (Q' = 0). 


passive system cannot have a negative imaginary 
part of to. The decay of the oscillation is caused 
by energy lost at the boundaries (recall that dp/dx 
= 0 does not represent complete reflection when 
the Mach number is nonnegligible). When the 
model was applied to small Mach number flows, 
the imaginary part of co approached zero for the 
passive case. 

Simple models of the fluctuating heat 
release rate were tested, in order to evaluate the 
model qualitatively. First, the heat release rate 
was described as a constant times the pressure 
perturbation, that is, Q' = c,p'. The specified 
value for the constant c, was large enough to 
cause a significant change in the eigenvalues of 
the system, relative to the passive case results. 
The new lowest frequency eigenvalue is shown 
as point C in Figure 4. The imaginary part of co 
decreased (i.e., the system became more 
unstable), when the fluctuating heat release rate 
was in phase with the pressure. This result is 
consistent with the Rayleigh criterion, and 
represents another important qualitative check of 
the model. The constant was then reduced to a 
third of its original value, in order to examine if 
the change in the growth rate of disturbances 
varied proportionally. The eigenvalue for the 
reduced constant is point D. The growth rate did 
not vary linearly with c,. 


Next, the fluctuating heat release rate was 
modeled as a purely imaginary constant times 
the pressure perturbation (i.e., the heat release 
rate was n/2 radians out of phase with the 
pressure). In this case two eigenvalues were 
obtained with frequencies similar to or less than 
those for the passive case. Both eigenvalues are 
labeled E in figure 4. This model for the 
fluctuating heat release rate should not add or 
subtract acoustic energy from the system, but 
only shift the phase. If there were no acoustic 
energy loss at the boundaries, we would expect 
the imaginary part of co to be zero for both the 
passive and active systems. Since there is a 
change in the acoustic energy lost at the 
boundaries as the frequency changes, there is a 
change in the decay rate of the disturbances 
when the heat release rate is modeled by this 
procedure. 

Point F represents the eigenvalue when the 
fluctuating heat release rate was in phase with 
the velocity perturbation, that is, Q' = c 2 u'. This 
representation of the heat release rate resulted in 
a smaller effect on the stability of the system than 
did the case Q' equal to a function of p'. This 
finding is to be expected, inasmuch as the 
standing velocity wave is far from being in phase 
with the standing pressure wave. If boundary 
conditions that caused the pressure and velocity 
to be in phase were selected, the results for this 
heat release rate model would parallel those for 
points C, D, and E. 

The last heat release rate model used the 
perturbations in chemistry described in the 
previous section. Point G is the eigenvalue for 
the first mode. Figure 6 shows the corresponding 
pressure and velocity standing waves. 

The second mode was also considered. 
Point B is the passive case (O' = 0) and point H 
the active case, that is, using the chemistry 
model. For the latter case the system was 
unstable (negative imaginary part of co); the 
pressure and velocity standing waves are given 
in figure 7. 

As final checks of the model, the effects of 
variable area were studied. The same flow 
conditions specified previously (see table 1 ) were 
used, but the duct area varied as follows: 

M=1 + 0.1x, (15) 

A(0) 
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Solid Line: normalized pressure 
Dashed Line: normalized velocity 


Figure 6. — Standing wave pattern for lowest frequency mode 
for constant area and active case 
(Q' from equation 12). 



Axial Position [m] 

Solid Line: normalized pressure 
Dashed Line: normalized velocity 

Figure 7. — Standing wave pattern for second lowest frequency 
mode for constant area and active case 
{Q' from equation 12). 


where A(0) is the combustor inlet area and A(x) 
the area at axial position x. Thus, there was a 
20% increase in area along the 2 meter long 
combustor. This change in area, of course, 
changes the mean flow quantities. Therefore, 
LSENS had to be used to calculate the new 
values for the mean flow variables and the partial 
derivatives 3Q/3p and 3Q/3T at each axial 
location. The eigenvalues for the first two modes 
are shown in figure 8. Points I and J represent 
the passive case (Q' = 0), while points K and L 


120 
100 

80 
60 

jf 40 
20 
0 

-20 

1000 1500 2000 2500 3000 

Re(w) 

Figure 8. — Eigenvalues for variable-area flow conditions and 
dynamic heat release rate models described in table 2. 
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represent the model of Q' using equation 12. 
Once again the amplification of the acoustic 
disturbances by the dynamic heat release rate 
caused the imaginary part of co to become more 
negative, indicating increasing system instability. 
For the variable-area case, the unsteady heat 
release due to the propane chemistry had a 
greater effect on the first mode than on the 
second. In contrast, the constant-area case 
displayed the opposite behavior (see figure 4). 

Conclusions 


A detailed chemical kinetics mechanism 
was successfully included in a linear stability 
model for one-dimensional reacting flow. By 
perturbing the chemical kinetic rates, the model 
was shown to predict qualitatively correct 
feedback between combustion and acoustics. 
When the pressure was in phase with the 
unsteady heat release rate, the system became 
less stable, which is consistent with the Rayleigh 
criterion. In addition, perturbations in chemical 
kinetic rates tended to make the system less 
stable. 

The model is being improved to incorporate 
fuel flow rate fluctuations, which must be included 
before the theoretical predictions can be 
compared with experimental results. Varying the 
fuel flow rate will produce perturbations in the 
local composition. By using detailed chemistry, 
the attendant perturbations in flame location and 
heat release rate can be modeled accurately. 
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